Tutorial Part IV — Getting started with animal movement

Introduction

In this tutorial, I will show you how to get started with (real) animal movement data. To put it simple, animal movement data are just points in time and space, i.e. re-locations of an individual can be provided as a data.frame with four columns. An identifier for the individual is only needed if you have different animals in one data set, i.e. if you have collared several individuals with the same collar or already appended several data sets:

Identifier	Timestamp	x_coordinate	y_coordinate
Fox_1	2022-01-01 08:10:04	13.47027	52.497802
Fox_1	2022-01-01 08:15:04	13.47015	52.497813

From Course2_SpatialR, you might remember that these can be imported as simple .txt or .csv files. Once you know in which coordinate reference system (CRS) your coordinates are stemming from, you can assign it to the data.frame, thereby you are creating a geo-referenced data set: a SpatialPointsDataframe - Object (if you work with the R-package {sp}) or an sf - Object (with package {sf}). Looking at these x- and y-coordinates, you might remember that they look like angular units, e.g. decimal degrees, and the x-coordinate refers to longitude and the y-coordinate to latitude. Hence, we are dealing with EPSG-code 4326. Once you have imported and georeferenced your data, you can plot and thoroughly check the data.

Checking your data is a mandatory prerequisite before any analysis (see Zuur et al. 2010). With your movement data, you might have outliers in space (e.g. when you test your collars in location A which might be hundreds of kilometers from your actual trapping site), your data might not have been regularly sampled because the animal was hiding and there was no access to the satellite (weak GPS-signal) etc etc. This uneven sampling can affect calculations of speed, for example, and therefore a series of packages have been developed to deal with those issues in movement data, like {move}.

The most prominent movement analyses comprise home range estimation, calculations of habitat preference, behavioral analyses (e.g. Hertel et al. 2021) and detection of movement syndromes (personality differences) (Michelangeli et al. 2021). Some of the packages listed in the next chapter are designed for these analyses.

TODO plot pic (fox_Q_20190113_sks.jpg “fox q”) funktioniert nicht

[TODO - exchange Q with Caro’s foxes, ref PhD C Scholz]

In the following, we will step by step start with loading and exploring movement data of female red fox (Vulpes vulpes) ‘Q von Stralau’, who had a small home range in Berlin, Germany. She was collared in January 2018 and had a very stable daily routine as shown by 4 (20) min relocation intervals. We will, however, only use the data of one month in this exercise, as data sets quickly get too big. Please note: the courses Course1_IntroR and Course2_SpatialR are obligatory for this tutorial.

Useful (web)sites and reading

For analysis of telemetry data: packages adehabitatHR, adehabitatLT, move,recurse, momentuHMM, moveHMM, ctmm, amt

Hertel, AG, Royauté, R, Zedrosser, A, Mueller, T. Biologging reveals individual variation in behavioural predictability in the wild. J Anim Ecol. 2021; 90: 723– 737. https://doi.org/10.1111/1365-2656.13406
Joo, R, Boone, ME, Clay, TA, Patrick, SC, Clusella-Trullas, S, Basille, M. Navigating through the r packages for movement. J Anim Ecol. 2020; 89: 248– 267. https://doi.org/10.1111/1365-2656.13116
Michelangeli, M., Payne, E., Spiegel, O., Sinn, D. L., Leu, S. T., Gardner, M. G., & Sih, A. (2021). Personality, spatiotemporal ecological variation and resident/explorer movement syndromes in the sleepy lizard. Journal of Animal Ecology, 00, 1– 14. https://doi.org/10.1111/1365-2656.13616
Scholz, C (2020). The ecology of red foxes (Vulpes vulpes) in anthropogenic environments. PhD thesis, FU Berlin.
Zuur, A.F., Ieno, E.N. and Elphick, C.S. (2010), A protocol for data exploration to avoid common statistical problems. Methods in Ecology and Evolution, 1: 3-14. https://doi.org/10.1111/j.2041-210X.2009.00001.x

Getting started

To follow the tutorial, you can either clone or download the repository at https://github.com/stephkramer/Course4_MoveQ or you create your own R-project, copy the raw data and type the code chunks into an R-Script. Please refer to the section on using R-projects in Course2_RSpatial.

If you start with your own R-project, I strongly recommend to use the d6-package Cedric Scherer provided. This package automatically sets up the ideal folder structure: https://github.com/EcoDynIZW/d6

In any case, my course folder has the following structure:

.
└── Course4_MoveQ
    ├── data-raw (contains the file with the GPS locations)
    ├── docs
    ├── output
    ├── plots
    └── R (contains the R-script)

Necessary packages to install and load

We first have to install the packages and load them before we can use the functions we need.

# package names:
pkgs = c("here", "sp", "sf", "dismo", "raster", "GISTools",  "rgdal", 
         "maptools", "rgeos","rgl","rasterVis", "adehabitatHR", "move", "tmap",
         "plotly", "circular","gganimate","moveVis", "ggmap", "maps", "mapproj",
         "viridis","dplyr", "devtools","lubridate", "patchwork",
         "colorspace","ragg", "ggtext","pdftools", "units",
         "glue", "cowplot", "tidyverse", "ggplot2", "Cairo") 

# install.packages(pkgs) # only run this line once for installing the packages!
# update.packages()

Tip of the day: If you have already installed some of the packages above, you can first check which ones are already installed, and save the ones not installed in an object called and only install the missing ones:

# my_packages <- pkgs[!(pkgs %in% installed.packages()[,"Package"])]
# my_packages
# if(length(my_packages) > 0) install.packages(my_packages)

Now load the packages:

library(here) #for easy directory management
library(sp)
library(dismo)
library(raster)
library(GISTools)
library(rgdal)    # retiring end of 2023 -> use stars/ terra /sf
library(maptools) # retiring end of 2023 -> use sp
library(rgeos)    # retiring end of 2023 -> use sf
library(rgl)
library(rasterVis)
library(adehabitatHR)  #adehabitatLT  #adehabitatMA                           
library(sf)
library(move)
library(plotly)
library(circular)
library(dplyr)
library(tmap)
library(viridis) 
library(gganimate)
library(ggplot2)
library(moveVis)
library(ggmap)
library(maps)  
library(mapproj)
library(devtools)
library(lubridate) 
library(here)
library(glue)
library(cowplot)
library(tidyverse)
library(Cairo)
library(colorspace)
library(ragg)
library(ggtext)
library(pdftools)
library(units)
library(ggplot2)
library(patchwork) # to combine plots

Set the working environment

Now that we are working inside an R-project, we can use the easy functionality of the {here} package, which is automatically pointing to your project’s root folder, e.g.:

here::here()

## [1] "C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub/Course4_MoveQ"

Hence, there is no need to use the function setwd() any more. Note: if it does not work, please close RStudio, go to your Explorer and double-click on the .Rproj file. Then, under ‘files’ (usually lower right panel) double-click on the R folder and open the script.

Load data

The movement data are stored in the data-raw subfolder. Let’s check which files are available.

lf <- list.files(path=here("data-raw"), full.names=TRUE) 
lf

## [1] "C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub/Course4_MoveQ/data-raw/geo-raw"        
## [2] "C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub/Course4_MoveQ/data-raw/tag5334_gps.txt"

The output lists two results, [1] is another subfolder, and [2] a text file. If you already know that your movement data contain e.g. ‘gps’ in their names or are stored as ‘.txt’ or ‘.csv’ files, you can directly search for those files with the pattern argument:

# check the difference, and note: full.names was set to FALSE
list.files(path=here("data-raw"),pattern='gps',full.names=FALSE)

## [1] "tag5334_gps.txt"

Now load the data file.

dat_anim <- read.table(file=lf[2], header=TRUE, fill=TRUE, sep=',')

Data check

Let’s have a look at the data. This is the typical data stored on e-obs collars:

dat_anim[1:5,] # recap: head(dat_anim) also works

##   key.bin.checksum tag.serial.number         start.timestamp longitude latitude
## 1       2135768252              5334 2018-10-30 04:00:00.000  13.47047 52.49499
## 2       2210577241              5334 2018-10-30 04:20:00.000  13.47157 52.49526
## 3       1824790132              5334 2018-10-30 04:40:00.000  13.47103 52.49525
## 4       4079427513              5334 2018-10-30 07:00:00.000  13.47109 52.49526
## 5        380027806              5334 2018-10-30 08:00:00.000  13.47105 52.49527
##   height.above.ellipsoid type.of.fix status used.time.to.get.fix
## 1                   76.3           3      A                   22
## 2                   81.6           3      A                   18
## 3                   66.5           3      A                   53
## 4                   73.5           3      A                   24
## 5                   71.1           3      A                   47
##          timestamp.of.fix battery.voltage fix.battery.voltage temperature
## 1 2018-10-30 04:00:23.000            3604                3389          10
## 2 2018-10-30 04:20:19.000            3616                3333          11
## 3 2018-10-30 04:40:54.000            3624                3374          13
## 4 2018-10-30 07:00:25.000            3659                3491          11
## 5 2018-10-30 08:00:48.000            3660                3470          13
##   speed.over.ground heading.degree speed.accuracy.estimate
## 1              1.64         233.31                    1.67
## 2              4.43         327.02                    3.79
## 3              0.07           0.00                    0.73
## 4              0.01           0.00                    0.65
## 5              0.04           0.00                    0.89
##   horizontal.accuracy.estimate
## 1                         7.94
## 2                        59.65
## 3                        15.36
## 4                         7.17
## 5                         8.45

For now, we will only work with few columns. tag.serial.number here refers to the individual identifier (collar ID). There are two columns with timestamps. start.timestamp is the preprogrammed time-interval. Then there is the timestamp.of.fix, which is the real time the GPS-location was recorded. This is usually a bit later, i.e.
( = start.timestamp + used.time.to.get.fix + 1 second), as it takes some time for the collar unit to connect to the satellite. The spatial info is stored in the columns longitude and latitude.

Before you can transform the data.frame ‘dat_anim’ into a georeferenced spatial object, you need to check whether there are missing locations in your data.frame, otherwise you will get an error message on transformation. This can happen if no GPS-signal could be recorded. Or - importantly - some collars are only activated when the animal is moving to save battery life. In that case, the missing GPS coordinates would correspond with the last position (= be the same). Depending on which analysis you want to do, e.g. define resting places, you might need to fill the missing positions again.

# if the latitude-entry is missing, the longitude value will also be missing
# so it is enough to only check the latitude
which(is.na(dat_anim$latitude)) #there are a lot of missing values in the locations

##   [1]   48   65   89   90   91  129  183  185  222  224  226  229  260  266  355
##  [16]  358  415  421  422  423  457  474  510  511  577  601  616  621  626  646
##  [31]  673  674  685  696  721  774  784  796  822  833  842  855  863  907  921
##  [46]  940  993 1006 1038 1070 1077 1078 1088 1126 1150 1159 1183 1188 1215 1238
##  [61] 1268 1272 1304 1330 1480 1516 1550 1552 1559 1600 1637 1680 1717 1721 1736
##  [76] 1738 1742 1772 1820 1824 1830 1867 1884 1973 1997 1999 2020 2023 2024 2025
##  [91] 2026 2064 2075 2076 2109 2162 2234 2313 2329 2374 2379 2396 2419 2453 2455
## [106] 2518 2551 2609 2614 2638 2687 2728 2745 2770 2776 2783 2800 2804 2805 2854
## [121] 2926 2927 2931 2935 2943 2950 2961 2963 2968 2990 2991 3043 3044 3052 3090
## [136] 3102 3114 3172 3225 3226 3305 3328 3361 3372 3373 3386 3409 3498 3516 3517
## [151] 3518 3520 3591 3622 3623 3624

dat_anim_na <- dat_anim[!is.na(dat_anim$latitude),] #alternatively, use complete.cases()

Make the crosscheck if there is missing data in longitude:

which(is.na(dat_anim_na$longitude)) # none

## integer(0)

Data preparation

We will now add some additional columns, where separate days (numbered from 1 to 365), the month (from 1 to 12) and the hour of the day are stored (1 - 23). This can be done with the package ´{lubridate}´. Check the vignette!

# define date-time format - the format is year-month-day_hour:min:sec:
dat_anim_na$start.timestamp <- ymd_hms(dat_anim_na$start.timestamp, tz="Europe/Berlin")

Now append the information:

dat_anim_na$yearday <- yday(dat_anim_na$start.timestamp)
dat_anim_na$month   <- month(dat_anim_na$start.timestamp)
dat_anim_na$hour    <- hour(dat_anim_na$start.timestamp)
dat_anim_na$kweek   <- week(dat_anim_na$start.timestamp)
# crosscheck with
# head(dat_anim_na)

In addition, we want to calulate the hours of sunset and sunrise as well as daylength. For this, we need to install a package that is still under development, i.e. which is not on CRAN. We therefore must download and install it locally:

# devtools::install_github("bgctw/solartime") # run only once
library(solartime)

For computing sunset and sunrise, the latitude/ longitude must be provided as well:

dat_anim_na$sunrise   <- computeSunriseHour(timestamp = dat_anim_na$start.timestamp,
                                            latDeg = dat_anim_na$latitude,
                                            longDeg = dat_anim_na$longitude)
dat_anim_na$sunset    <- computeSunsetHour(dat_anim_na$start.timestamp,
                                           dat_anim_na$latitude,
                                           dat_anim_na$longitude)
dat_anim_na$daylength <- computeDayLength(dat_anim_na$start.timestamp,dat_anim_na$latitude)
dat_anim_na$daytime   <- computeIsDayByLocation(dat_anim_na$start.timestamp,
                                                dat_anim_na$latitude,
                                                dat_anim_na$longitude)
# head(dat_anim_na)

Data transformation

Now we can transform the data into sf and sp objects and assign the coordinate reference system:

# transform into spatial simple feature sf object
mydf_sf <- st_as_sf(x = data.frame(dat_anim_na),
                       coords = c("longitude", "latitude"),
                       crs = 4326,
                       sf_column_name = "geometry" )

# transform into SpatialPointsDataFrame  - for crosschecking
mydf_sp <- as(mydf_sf, "Spatial")

And then we project the reference system from angular units to a planar coordinate reference system in meters:

# transform CRS to projected one in meter distance units
mydf_sf_trans <-  st_transform(mydf_sf, 5631 )  # EPSG-code Pulkovo 
mydf_sp_trans <-  spTransform(mydf_sp, CRS("+init=epsg:5631"))

Data exploration

Plot the data

This is recap from Course2. A quick plot from the simple data.frame

# latitude is y-Axis, longitude = x-axis in cartesian coordinate system
plot(dat_anim_na$latitude ~ dat_anim_na$longitude,
       pch='.',
       col=as.numeric(dat_anim_na$daytime)+1 # + 1 because color 0 is transparent
     )

A quick plot of the movement path, using only the first 50 rows

plot(dat_anim_na$latitude[1:50] ~ dat_anim_na$longitude[1:50], type= 'l')

Much nicer with ggplot and the projected sf object:

baseplot <- ggplot(data = mydf_sf_trans) +
              geom_sf(aes(color = daytime),size=0.01, alpha=0.5)   +  
              xlab("Longitude") + ylab("Latitude") +
              ggtitle("Telemetry data")
baseplot

We can even add contour lines of the point density: TODO - error here - duplicated aes

densplot <- baseplot + stat_density_2d(
                         aes(x = st_coordinates(mydf_sf_trans)[,1], 
                             y = st_coordinates(mydf_sf_trans)[,2],
                         alpha = ..level..,
                         color = after_scale(colorspace::darken(color, .5, space = "HLS")),
                         fill  = after_scale(colorspace::lighten(color, .4, space = "HLS"))), 
                         geom  = "polygon",
                         color = "white" ,
                         size  = .02,
                         n     = 500,
                         bins  = 5,
                         adjust = 3,
                         data  = mydf_sf_trans)  
densplot

# example of how to save the plot to your file
# ggsave(here("plots","plot_pointdensity.png"), dpi = 800, height = 6, width = 6)

We can also add background using ggmap. Check here for a quick start to ggmap: https://www.nceas.ucsb.edu/sites/default/files/2020-04/ggmapCheatsheet.pdf

# define abounding box for the plot
b <- c(13.460,52.485,13.490,52.500) # make large enough to plot all data!

t <- get_map(location = b, maptype = "terrain", source='osm', zoom = 15)

(ggmap(t) 
  + geom_point(aes(x = longitude, y = latitude), 
               data = data.frame(dat_anim_na), 
               colour = as.numeric(dat_anim_na$daytime)+5, 
               size = 0.05)
  + labs(x="Latitude", y="Longitude", title="telemetry locations") 
)

Last but not least and much quicker, the ´{tmap}´package can help, but you need an sf/sp-object:

tmap_mode(mode = "view")
  tm_shape(shp = mydf_sf_trans)  + tm_dots(size = 0.01, 
                                           col = 'daytime',  
                                           alpha=0.5)

Have a deep look at the data: did the fox really swim? (points in Lake Rummelsburg). Or are these outliers? No: the lake was frozen! Note: you really need to know your animal and the area!

END

Adding the session info can be very helpful when going back to old scripts or using scripts of others:

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: i386-w64-mingw32/i386 (32-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_Germany.1252  LC_CTYPE=English_Germany.1252   
## [3] LC_MONETARY=English_Germany.1252 LC_NUMERIC=C                    
## [5] LC_TIME=C                       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] solartime_0.0.3     patchwork_1.1.1     units_0.7-2        
##  [4] pdftools_3.0.1      ggtext_0.1.1        ragg_1.2.1         
##  [7] colorspace_2.0-2    Cairo_1.5-14        forcats_0.5.1      
## [10] stringr_1.4.0       purrr_0.3.4         readr_2.1.1        
## [13] tidyr_1.1.4         tibble_3.1.6        tidyverse_1.3.1    
## [16] cowplot_1.1.1       glue_1.5.1          lubridate_1.8.0    
## [19] devtools_2.4.3      usethis_2.1.5       mapproj_1.2.7      
## [22] maps_3.4.0          ggmap_3.0.0         moveVis_0.10.5     
## [25] gganimate_1.0.7     viridis_0.6.2       viridisLite_0.4.0  
## [28] tmap_3.3-2          dplyr_1.0.7         circular_0.4-93    
## [31] plotly_4.10.0       ggplot2_3.3.5       move_4.1.6         
## [34] geosphere_1.5-14    sf_1.0-4            adehabitatHR_0.4.19
## [37] adehabitatLT_0.3.25 CircStats_0.2-6     boot_1.3-28        
## [40] adehabitatMA_0.3.14 ade4_1.7-18         deldir_1.0-6       
## [43] rasterVis_0.51.1    lattice_0.20-45     rgl_0.108.3        
## [46] rgdal_1.5-27        GISTools_0.7-4      rgeos_0.5-9        
## [49] MASS_7.3-54         RColorBrewer_1.1-2  maptools_1.1-2     
## [52] dismo_1.3-5         raster_3.5-11       sp_1.4-6           
## [55] here_1.0.1         
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2              tidyselect_1.1.1        htmlwidgets_1.5.4      
##   [4] grid_4.1.2              munsell_0.5.0           codetools_0.2-18       
##   [7] av_0.6.0                gifski_1.4.3-1          withr_2.4.3            
##  [10] highr_0.9               knitr_1.36              rstudioapi_0.13        
##  [13] wk_0.5.0                labeling_0.4.2          RgoogleMaps_1.4.5.3    
##  [16] farver_2.1.0            rprojroot_2.0.2         vctrs_0.3.8            
##  [19] generics_0.1.1          xfun_0.28               qpdf_1.1               
##  [22] R6_2.5.1                slippymath_0.3.1        rmdformats_1.0.3       
##  [25] isoband_0.2.5           bitops_1.0-7            cachem_1.0.6           
##  [28] assertthat_0.2.1        scales_1.1.1            gtable_0.3.0           
##  [31] lwgeom_0.2-8            processx_3.5.2          rlang_0.4.12           
##  [34] systemfonts_1.0.3       lazyeval_0.2.2          dichromat_2.0-0        
##  [37] hexbin_1.28.2           broom_0.7.10            s2_1.0.7               
##  [40] yaml_2.2.1              abind_1.4-5             modelr_0.1.8           
##  [43] crosstalk_1.2.0         backports_1.4.1         gridtext_0.1.4         
##  [46] tools_4.1.2             bookdown_0.24           ellipsis_0.3.2         
##  [49] jquerylib_0.1.4         proxy_0.4-26            sessioninfo_1.2.2      
##  [52] Rcpp_1.0.7              plyr_1.8.6              base64enc_0.1-3        
##  [55] progress_1.2.2          classInt_0.4-3          ps_1.6.0               
##  [58] prettyunits_1.1.1       pbapply_1.5-0           tmaptools_3.1-1        
##  [61] zoo_1.8-9               haven_2.4.3             fs_1.5.2               
##  [64] leafem_0.1.6            magrittr_2.0.1          data.table_1.14.2      
##  [67] magick_2.7.3            reprex_2.0.1            mvtnorm_1.1-3          
##  [70] pkgload_1.2.4           hms_1.1.1               evaluate_0.14          
##  [73] XML_3.99-0.8            leaflet_2.0.4.1         jpeg_0.1-9             
##  [76] readxl_1.3.1            gridExtra_2.3           testthat_3.1.1         
##  [79] compiler_4.1.2          KernSmooth_2.23-20      crayon_1.4.2           
##  [82] htmltools_0.5.2         tzdb_0.2.0              DBI_1.1.1              
##  [85] tweenr_1.0.2            dbplyr_2.1.1            cli_3.1.0              
##  [88] parallel_4.1.2          pkgconfig_2.0.3         foreign_0.8-81         
##  [91] terra_1.4-22            xml2_1.3.3              bslib_0.3.1            
##  [94] rvest_1.0.2             callr_3.7.0             digest_0.6.29          
##  [97] rmarkdown_2.11          cellranger_1.1.0        leafsync_0.1.0         
## [100] curl_4.3.2              rjson_0.2.20            lifecycle_1.0.1        
## [103] jsonlite_1.7.2          desc_1.4.0              askpass_1.1            
## [106] fansi_0.5.0             pillar_1.6.4            fastmap_1.1.0          
## [109] httr_1.4.2              pkgbuild_1.3.0          remotes_2.4.2          
## [112] png_0.1-7               leaflet.providers_1.9.0 class_7.3-19           
## [115] stringi_1.7.6           sass_0.4.0              textshaping_0.3.6      
## [118] latticeExtra_0.6-29     stars_0.5-5             memoise_2.0.1          
## [121] e1071_1.7-9